1 Introduction

Lab book for analyses using hierachal computational modelling to identify paramters that define the best model of learning as it applies to fear conditioning acquisition and extinction using FLARe fear conditioning data. Long abstract, justification and analysis plan found in prelim manuscript here

In short:

1.1 Aims

Identify model of learning based on a priori hypotheses that best fits the trajectories of fear relevant learning in our FLARe dataset
- Use all first week data from Validation, app TRT, lab TRT, Pilot, Headphones (n = 223 after exclusions)
- Include Acquisition, extinction (trajectories representing fear learning and treatment)
- Identify parameters that define these trajectories
  - e.g. Learnign rate, plateau, first ambiguous trial etc.
Cross validate best fitting model in TEDS data
Are these parameters associated with other emasures of indsividual differences in our datasets?
- Personality (Neuroticism)
- Current anxiety symptoms (GAD-7) - equivalent of baseline symptoms (Chris + Meg analyses)
- Lifetime / trait anxiety (STAI / ASI - FLARe analyses)
- Current depression symptoms (PHQ-9) - equivalent of baseline symptoms (Chris + Meg analyses)
- Interpretation biases (IUS, ASSIQ - FLARe analyses)
- SES (Meg IAPT: benefits, employment etc)
- Gender (Meg analyses)
- Emotion regulation profile (potentially LCA based?)

1.2 Impact and relevance

Evidence from both human (Richter et al., 2012) and rodent (Galatzer-Levy, Bonanno, Bush, & LeDoux, 2013) studies suggest that trajectories of how we learn and extinguish fear differ between individuals. Different trajectories of fear and extinction have also been found using fear conditioning studies (e.g. Duits et al., 2016), a good model for the learning of, and treatment for, fear and anxiety disorders. It is likely that these trajectories of fear extinction might predict outcomes in exposure-based cognitive behavioural therapy (Kindt, 2014). 
 
Identifying parameters that predict individual trajectories of fear learning and extinction will enable us to harness fear conditioning data more effectively to aid in understanding mechanisms underlying the development of and treatment for anxiety disorders. With more accurate models of these processes, the potential to use fear conditioning paradigms to predict those most at risk of developing an anxiety disorder, and those who might respond best to exposure-based treatments, greatly improves.

1.3 Useful references

Sutton and Barto Reinforcement Learning - Textbook on reinforcement learning
Anxiety promotes memory for mood-congruent faces but does not alter loss aversion (Charpentier…Robinson, 2015) - Good example of a sensitivity learning parameter
Hypotheses About the Relationship of Cognition With Psychopathology Should be Tested by Embedding Them Into Empirical Priors (Moutoussist et al., 2018) - Including variables of interest (e.g. anxiety) in the model

Toby Wise has just submitted an aversive learning paper incorporating beta probability distributions in the best model for uncertain learning paramters etc.

A copy of this is

1.4 Analysis plan

Define set of a priori models moving from simple to more complex
- Some paramters to include:
  - Rate of learning (sometimes with punishment reinforcement)
  - Sensitivity to punishment
  - Pre-existing anxiety
  - SES? Gender?
Run each model and compare fit in FLARe pre TEDS data
- Use Log likelihood and BIC etc.
Select best fitting model
Extract individual data for learning parameters from this model and see what factors best predict it
- Anxiety (if anxiety isnt best as part of the model)
- Interpretation biases
- Tolerance of uncertanty
- Cognitive emotional control
- emotional attentional control
- SES?
- Gender?
Run all models again in FLARe TEDS
- Decide if the same model best fits the data again.
- See if we get similar results from the parameter prediction

Will use a combination of R.Version(3.5.1), RStan (Version 2.18.2, GitRev: 2e1f913d3ca3) and hBayesDM package in R (3.5.1) Ahn, W.-Y., Haines, N., & Zhang, L. (2017). Revealing neuro-computational mechanisms of reinforcement learning and decision-making with the hBayesDM package. Computational Psychiatry, 1, 24-57., which uses RStan

1.5 Modelling notes

1.5.1 Intuition

Discussion with Vince Valton and Alex Pike about the best way to fit this model. As the observed outcomes (expectancy ratings) are non binary and are related to eachother (i.e. as you become more likely to select 9, you become less likely to select 1) we should consider each trial for each person for each stimulus as a constantly updating beta distribution. so you might see a pattern like this for the CS+ in acq for example.

So, best model is likely to be one using beta distributions that show the probability distribution for each rating.

We can use sufficient parameters to describe these (i.e. mean / sd or possibly the mode)

scaling

We can scale the beta by how aversive participants find the shock. i.e. it might update their learning as if there was .5 a shock or 1.5 of a shock depending on their own sensitivity to the aversiveness / punishment.

alpha

generalisation

We can do this with a single beta distribution for each phase (collapsing over the two stimuli). This would be akin to a per phase generalisation paramaterer in that it will be smaller if they tend to choose the same expectancy for both stimuli and larger if they tend to choose very differently for both stimuli.

However, because these variables are not really equivalent (i.e the reinforcement rate is different for both, and we use this in the model)

So instead we can create a paramater which is the value of cs- weighted by some value of the cs+. How much each individual weights by the Cs+ can be freely estimated by the model and can be the generalisation paramter.

So this would be vminus = vminus + (w)vplus (where the w paramter is the freely estimated paramter per person)

per stimulus We probably want to model cs+ and cs- separately too - so have a beta distribution characterised by sufficient parameters for each.

per trial

All of the above can then also be done with updating per trial.

leaky beta

we also need a model that incorporates ‘leak’. i.e. learning leak - likely that participants will update more based on more recent trials and learn less from the more distant trials as time progresses. See Toby’s paper for more.

uncertainty

We should consider incorportating a paramter that maps to participant uncertainty about outcomes.

anxiety

Might be worth incorporating this as a model paramater / feature. Read this for more.

Hypotheses About the Relationship of Cognition With Psychopathology Should be Tested by Embedding Them Into Empirical Priors (Moutoussist et al., 2018)

1.5.2 Terminology

V == ‘value’. Baasically a paramter that is about the salience of the stimulus at any given point.
alpha == ‘learning rate’. A parameter that describes how sensitive people are to updating their learning. So a fast learning rate means that learning on any given trial is weighted more based on the trials immediatly preceding than past ones, and a slow learning rate means that all past events influence learning more evenly. Alex’s tennis analogy is good here (Federer - stable player, can predict a win based on all matches; Murray - volatile player; his last match is best predictor of next match performance).

1.5.3 Example of beta distributions

and how they change depending on whether you change the beta or alpha paramters.

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

##  [1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65
## [15] 0.70 0.75 0.80 0.85 0.90 0.95 1.00

## [1] "stable beta, increasing alpha"

## [1] "stable alpha, increasing beta"

1.5.4 Models to write / run

Will probably do all per trial. Will do an early sensitivity check to confirm this.

Single beta, no scaling
Single beta, no scaling per trial.
*** At this point, compare the two above. Ensure the per trial fits better, and if it does then do all below per trial***
“” scaled
Single beta Single alpha reinforcement learning model (estimate both the beta and the alpha i.e. learning rate)
Single beta single alpha reinforcement learning with mean + sd for the beta estimate as a paramter
Beta per stimulus
Beta per stimulus + generalisation paramter (Vminus = vminus + wvplus)
Leaky beta
Leaky beta + uncertainty
Leaky beta + uncertainty + anxiety

2 Analyses

2.1 Preliminary

2.1.1 Compare a priori to data

2.1.1.1 Simulate different learning rates

only doing this ‘accurately’ for the acquisition CS+, as the simulations require probability. I am using contingency for this (0.75). If set for 0 for all other phases and stimuli then it looks as if the learning should be flat regardless of alpha. We expect in reality that this probability will vary between people and will be unlikely to be zero. So test 12 and 18 trials with a probability of 0.5 and 0.2 as well.

## [1] "Simulated learning rates. 12 trials; probability = 0.75 (CSp acq contingency) \n"

## [1] "Simulated learning rates. 12 trials; Probability = 0.5\n"

## [1] "Simulated learning rates. 12 trials; Probability = 0.2\n"

## [1] "Simulated learning rates. 18 trials; Probability = 0.5\n"

## [1] "Simulated learning rates. 18 trials; Probability = 0.2\n"

2.1.1.2 Plot subset of trajectories in flare

## Warning: package 'data.table' was built under R version 3.5.2

## 
## Attaching package: 'reshape2'

## The following objects are masked from 'package:data.table':
## 
##     dcast, melt

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:data.table':
## 
##     between, first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

2.1.2 Set up

These use Alex Pikes RStan script with minor modification to make it punishment only to see if it runs. Testing that the approach works with the current data set up etc.

The settings for the script are below, including stan chain paramters and directory set up.

This loads the libraries and source files needed to run this script, and sets up RStan

## Warning: package 'StanHeaders' was built under R version 3.5.2

2.1.3 Try RStan

See if the basic punishment only learning model for the CS+ and CS- works with the FLARe master data

2.1.3.1 Run the 8schools check

From the rstan github

This is to check that all is compiling and working and to give and idea of data format etc.

2.1.3.2 Adjust dataframe

load in the week 1 app and lab data for FLARe pilot, TRT and headphones studies. Make it long form.

Try with acquisition data first. This is formatted with no column names, with no missing data.

Derive the n parameter for both files and check these match

stanname='punish_only.stan'
minus_name <- 'bayes_acq_minus.csv'
plus_name <- "bayes_acq_plus.csv"

stanfile <- file.path(scriptdir, stanname)
minusfile <- file.path(datadir,minus_name)
plusfile <- file.path(datadir,plus_name)


minus <- fread(minusfile,data.table=F)
plus <-fread(plusfile,data.table=F)

nacqm <- dim(minus)[1]
nacqp <- dim(plus)[1]


## check that these match and create nsub variable for RStan

if (nacqm == nacqp) {
  print('subject number match')
  nsub <- nacqm
  
  print(paste('nsub set to',nsub,sep=" "))
} else {
  print('WARNING: subject number does not match. Check master dataset')
}

## [1] "subject number match"
## [1] "nsub set to 342"

# check the file format is ok

minus[1:2,]

plus[1:2,]

# create the n trials variable for RStan

ntrials <- 12

The expectancy rating datasets look like they are formatted fine and ntrials and nsub variables should exist.

2.1.3.2.1 Create scream data

Need to go back to stage zero and keep scream yes/no as a variable. For now to see if this runs create simulated version for the CS+. CS- will remain the same.

screamMinus <- matrix(0L,nrow=nsub, ncol=ntrials)

# Initialise plus dataset in the same way, but make the first trial 1 for everyone, then add 8 additional random 1's per person. Do this in four random patterns to mimic the real data

sc1 <- c(1,1,0,1,0,0,1,1,1,1,1)
sc2 <- c(0,1,1,1,0,0,1,1,1,1,1)
sc3 <- c(1,1,1,0,1,0,1,0,1,1,1)
sc4 <- c(1,0,1,1,0,0,1,1,1,1,1)


screamPlus <- matrix(0L,nrow=nsub, ncol=ntrials)

screamPlus[,1] <- 1


# for (n in 1:dim(screamPlus)[1]) {
#   print(n)
#   screamPlus[n,2:12] <- sample(patts,1,replace=T)
# }


for (n in 1:dim(screamPlus)[1]) {
  
  a <- sample(c(1,4),1)
  
  if (a == 1) {
    screamPlus[n,2:12] <- sc1
  } else if (a == 2) {
    screamPlus[n,2:12] <- sc2
  } else if (a == 3){
    screamPlus[n,2:12] <- sc3
  } else {
    screamPlus[n,2:12] <- sc4
  }
}

2.1.3.2.2 make rating data binary

for now to see if stan runs using bernoulli-logit function make binary resposnes from expectancy i.e. >=4.5 ==1, <= 4.5 ==0.

binarise <- function(x) {
  ifelse(x >= 4.5,1,0)
}

minusb <- data.frame(apply(minus,2,function(x) binarise(x)))

plusb <- data.frame(apply(plus,2,function(x) binarise(x)))

2.1.3.3 Set up procedure to create and sync models.

This directs to my local machine here /Users/kirstin/Dropbox/SGDP/FLARe/FLARe_MASTER/Projects/Hierachal_modelling/Scripts and is remotely linked to the github repository here.

2.1.3.3.1 Make sure the most up to date stan file is in the remote repo


git pull Bayes_modelling

## Already up to date.

2.1.3.3.2 check existing paramters

Unhash this if you want to check what the model looks like within the notebook.

stanname="punish_only.stan"
scriptdir="/Users/kirstin/Dropbox/SGDP/FLARe/FLARe_MASTER/Projects/Hierachal_modelling/Scripts"

#cat $scriptdir/$stanname

2.1.3.3.3 Make any changes

use echo to push these to the new file if you want to make changes from here.


## initialise bash directory and filename

stanname="punish_only.stan"
scriptdir="/Users/kirstin/Dropbox/SGDP/FLARe/FLARe_MASTER/Projects/Hierachal_modelling/Scripts"


#echo "<any changes here>" > $scriptdir/$stanname

2.1.3.4 Now try run stan

unhash this to run experimental script that checked if stan runs. This was mostly to check data formatting and installation / compilation etc.

flare_data<-list(ntrials=ntrials,nsub=nsub,includeTrial = rep(1,ntrials), screamPlus=t(screamPlus),screamMinus=t(screamMinus),
                 ratingPlus=t(plusb),ratingMinus=t(minusb))

#flare_fit <- stan(file = stanfile, data = flare_data, iter=chain_iter, chains = chain_n) #add working dir?

#save(flare_fit, file=file.path(datadir,'flare_fit_test'))

#traceplot(flare_fit,'lp__')

# extract fit data
#summary_flare<- summary(flare_fit)

# extract model summary data

#flare_loglike<- extract_log_lik(flare_fit, parameter_name = "loglik", merge_chains = TRUE)

view the fit information

#summary_flare

extract the loglikelihood using loo

#loo(flare_loglike)

so, good news is this all works. So preliminary check a success. Next need to consider the appropriate model.

2.2 Create stan friendly datasets

2.2.1 notes

We need to rescale our dataset here to be between 0 and 1.

Importantly, because we are using the proportion of trials that are not reinforced as a known paramter for statistical reasons (we don’t want a proportion of .75 and 1, better to have .25 and 0), we have made our rescaled expectancy values as 1 - rescaled(x). This means that we will still be able to interpret the results in the expected way (i.e. higher rating is greater expectation of the outcome).

2.2.2 rescale data

rescale the 1-9 expectancy values to be on a 0-1 scale.

stan cannot deal with the extreme limit of the beta, so make the rescaled limits just above 0 and below one

library(scales)

# rescale and flip so that we are effectively rating the expectation that they WILL NOT hear a scream to match stan

minus_scaled <- data.frame(apply(minus,2,function(x) 1-rescale(x, to=c(0.00001,0.99999))))
plus_scaled <- data.frame(apply(plus,2,function(x) 1-rescale(x, to=c(0.00001,0.99999))))

2.2.3 create proportion screams data

This is a vector containing the absolute number of trials where no scream occurred for each stimulus. As there was a 75% reinforcement rate for the CS+ (9/12 trials), this is a vector of ’3’s. For the CS-, no trials were reinforced so is a vector of ’12’s

No_scream_p <- rep(3,nsub)
No_scream_m <- rep(12,nsub)

2.2.3.1 Create scream data

Create datasets for the acquisition CS- and extinction CS+ and CS- reflecting that no screams occurred at all. Then use the pattern id variable to create a dataset for the acquisition CS+ indicating when a scream occurred for each participant.

## Create the no scream daatsets for all

screamMinus <- matrix(0L,nrow=nsub, ncol=ntrials)

# Initialise plus dataset in the same way, but make the first trial 1 for everyone, then add 8 additional random 1's per person. Do this in four random patterns to mimic the real data

sc1 <- c(1,1,0,1,0,0,1,1,1,1,1)
sc2 <- c(0,1,1,1,0,0,1,1,1,1,1)
sc3 <- c(1,1,1,0,1,0,1,0,1,1,1)
sc4 <- c(1,0,1,1,0,0,1,1,1,1,1)


screamPlus <- matrix(0L,nrow=nsub, ncol=ntrials)

screamPlus[,1] <- 1


# for (n in 1:dim(screamPlus)[1]) {
#   print(n)
#   screamPlus[n,2:12] <- sample(patts,1,replace=T)
# }


for (n in 1:dim(screamPlus)[1]) {
  
  a <- sample(c(1,4),1)
  
  if (a == 1) {
    screamPlus[n,2:12] <- sc1
  } else if (a == 2) {
    screamPlus[n,2:12] <- sc2
  } else if (a == 3){
    screamPlus[n,2:12] <- sc3
  } else {
    screamPlus[n,2:12] <- sc4
  }
}

2.3 Model 1: single beta no scaling

2.3.1 notes

Because we use the 1-rescaled expectancy data, no need to try and invert to reinforcement paramters here. As a result we need the stan model to simply be:

alphaPlus[p] =  nothingPlus[p]/ntrials;
alphaMinus[p] =  nothingMinus[p]/ntrials;

2.3.2 run Alex Pike’s stan script for non scaled beta model.

here we try to estimate the alpha paramter of the beta distribution per trial per person per stimulus. (i.e. you have two sufficient paramters for each beta dist, the alpha and beta. we want to estimate the alpha - ).

Eventually we will scale these by the actual ‘value’ of the scream for each person per trial.

Using data loaded in from preliminary tests above.

so this is a beta value per person (assuming the underlying process for the plus and minus are the same)

stanname='beta_noscaling.stan'

stanfile <- file.path(scriptdir, stanname)

flare_data<-list(ntrials=ntrials,nsub=nsub,nothingPlus = No_scream_p, nothingMinus=No_scream_m,ratingsPlus=plus_scaled,ratingsMinus=minus_scaled)

flare_fit <- stan(file = stanfile, data = flare_data, iter=chain_iter, chains = chain_n) #add working dir?

## 
## SAMPLING FOR MODEL 'beta_noscaling' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001121 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.21 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 12.2827 seconds (Warm-up)
## Chain 1:                4.55183 seconds (Sampling)
## Chain 1:                16.8345 seconds (Total)
## Chain 1:

save(flare_fit, file=file.path(datadir,'flare_fit_test'))

traceplot(flare_fit,'lp__')

# extract fit data
summary_flare<- summary(flare_fit)

# extract model summary data

#flare_loglike<- extract_log_lik(flare_fit, parameter_name = "loglik", merge_chains = TRUE)

2.4 Model 3: single beta scaled

2.4.1 notes

Simple alteration of the first model. We estimate a scaling parameter per person over all trials and apply this to alpha component per participant.

2.4.2 run Alex Pike’s stan script for scaled beta model.

Eventually we will scale these by the actual ‘value’ of the scream for each person per trial.

Using data loaded in from preliminary tests above.

so this is a beta value per person (assuming the underlying process for the plus and minus are the same)

stanname='beta_scaling.stan'

stanfile <- file.path(scriptdir, stanname)

flare_data<-list(ntrials=ntrials,nsub=nsub,nothingPlus = No_scream_p, nothingMinus=No_scream_m,ratingsPlus=plus_scaled,ratingsMinus=minus_scaled)

flare_fit <- stan(file = stanfile, data = flare_data, iter=chain_iter, chains = chain_n) #add working dir?

## 
## SAMPLING FOR MODEL 'beta_scaling' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001412 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 10.6291 seconds (Warm-up)
## Chain 1:                9.68512 seconds (Sampling)
## Chain 1:                20.3142 seconds (Total)
## Chain 1:

save(flare_fit, file=file.path(datadir,'flare_fit_test'))

traceplot(flare_fit,'lp__')

# extract fit data
summary_flare<- summary(flare_fit)

# extract model summary data

#flare_loglike<- extract_log_lik(flare_fit, parameter_name = "loglik", merge_chains = TRUE)

2.5 Model 4: reinforcement learning model, single beta, single alpha

2.5.1 notes

this model includes an alpha learning paramater per person estimating their learning rate and updating based on it. This model needs a dataset that indicates whether a scream occurred for each trial instead of the proportion of times no scream occurred.

2.5.2 run Alex Pike’s stan script for scaled beta model.

Eventually we will scale these by the actual ‘value’ of the scream for each person per trial.

Using data loaded in from preliminary tests above.

so this is a beta value per person (assuming the underlying process for the plus and minus are the same)

stanname='beta_withRL.stan'

stanfile <- file.path(scriptdir, stanname)

flare_data<-list(ntrials=ntrials,nsub=nsub,screamPlus = t(screamPlus), screamMinus= t(screamMinus),ratingsPlus=t(plus_scaled),ratingsMinus=t(minus_scaled))

flare_fit <- stan(file = stanfile, data = flare_data, iter=chain_iter, chains = chain_n) #add working dir?

## 
## SAMPLING FOR MODEL 'beta_withRL' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.005059 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 50.59 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 49.5232 seconds (Warm-up)
## Chain 1:                32.1446 seconds (Sampling)
## Chain 1:                81.6678 seconds (Total)
## Chain 1:

save(flare_fit, file=file.path(datadir,'flare_fit_test'))

traceplot(flare_fit,'lp__')

# extract fit data
summary_flare <- summary(flare_fit)

# extract model summary data

#flare_loglike<- extract_log_lik(flare_fit, parameter_name = "loglik", merge_chains = TRUE)

2.6 Model 5: reinforcement learning model, single beta, single alpha with mean SD defining beta shape

2.6.1 notes

2.6.2 run stan script

stanname='beta_meansd_RL.stan'

stanfile <- file.path(scriptdir, stanname)

flare_data<-list(ntrials=ntrials,nsub=nsub,screamPlus = t(screamPlus), screamMinus=t(screamMinus),ratingsPlus=t(plus_scaled),ratingsMinus=t(minus_scaled))

flare_fit <- stan(file = stanfile, data = flare_data, iter=chain_iter, chains = chain_n) #add working dir?

## 
## SAMPLING FOR MODEL 'beta_meansd_RL' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.009465 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 94.65 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 42.9256 seconds (Warm-up)
## Chain 1:                45.4189 seconds (Sampling)
## Chain 1:                88.3446 seconds (Total)
## Chain 1:

## Warning: There were 33 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

save(flare_fit, file=file.path(datadir,'flare_fit_test'))

traceplot(flare_fit,'lp__')

# extract fit data
summary_flare<- summary(flare_fit)

# extract model summary data

#flare_loglike<- extract_log_lik(flare_fit, parameter_name = "loglik", merge_chains = TRUE)

3 questions

why does the beta RL mean sd not run unless beta is bounded between 0 and 0.0001??

4 to do

add a beta for each stimulus
a version with mode instead of mean (check for equations?)
per trial for beta_noscaling

5 Push any updates to github

5.0.0.0.1 Any push the updates to github

Unhash the series below if you made any changes.


## initialise bash directory and filename
stanname="punish_only.stan"
scriptdir="/Users/kirstin/Dropbox/SGDP/FLARe/FLARe_MASTER/Projects/Hierachal_modelling/Scripts"


## stage
#git add $scriptdir/$stanname


## push 

#git push Bayes_modelling

Hierachal computational modelling - FLARe

1 Introduction

1.1 Aims

1.2 Impact and relevance

1.3 Useful references

1.4 Analysis plan

1.5 Modelling notes

1.5.1 Intuition

1.5.2 Terminology

1.5.3 Example of beta distributions

1.5.4 Models to write / run

2 Analyses

2.1 Preliminary

2.1.1 Compare a priori to data

2.1.1.1 Simulate different learning rates

2.1.1.2 Plot subset of trajectories in flare

2.1.2 Set up

2.1.3 Try RStan

2.1.3.1 Run the 8schools check

2.1.3.2 Adjust dataframe

2.1.3.2.1 Create scream data

2.1.3.2.2 make rating data binary

2.1.3.3 Set up procedure to create and sync models.

2.1.3.3.1 Make sure the most up to date stan file is in the remote repo

2.1.3.3.2 check existing paramters

2.1.3.3.3 Make any changes

2.1.3.4 Now try run stan

2.2 Create stan friendly datasets

2.2.1 notes

2.2.2 rescale data

2.2.3 create proportion screams data

2.2.3.1 Create scream data

2.3 Model 1: single beta no scaling

2.3.1 notes

2.3.2 run Alex Pike’s stan script for non scaled beta model.

2.4 Model 3: single beta scaled

2.4.1 notes

2.4.2 run Alex Pike’s stan script for scaled beta model.

2.5 Model 4: reinforcement learning model, single beta, single alpha

2.5.1 notes

2.5.2 run Alex Pike’s stan script for scaled beta model.

2.6 Model 5: reinforcement learning model, single beta, single alpha with mean SD defining beta shape

2.6.1 notes

2.6.2 run stan script

3 questions

4 to do

5 Push any updates to github

5.0.0.0.1 Any push the updates to github